function plotCubic(obj,funcM)
f = figure;
subplot(2,2,1)
contourf(obj.unitGrid.posList,obj.unitGrid.posList,...
    reshape(funcM(:,:,1),length(obj.unitGrid.posList),length(obj.unitGrid.posList)))
title('z=0');

subplot(2,2,2)
contourf(obj.unitGrid.posList,obj.unitGrid.posList,...
    reshape(funcM(:,:,ceil(length(obj.unitGrid.posList)/3)),length(obj.unitGrid.posList),length(obj.unitGrid.posList)))
title('z=L/3');

subplot(2,2,3)
contourf(obj.unitGrid.posList,obj.unitGrid.posList,...
    reshape(funcM(:,:,ceil(2*length(obj.unitGrid.posList)/3)),length(obj.unitGrid.posList),length(obj.unitGrid.posList)))
title('z=2L/3');

subplot(2,2,4)
contourf(obj.unitGrid.posList,obj.unitGrid.posList,...
    reshape(funcM(:,:,end),length(obj.unitGrid.posList),length(obj.unitGrid.posList)))           
title('z=L');

colormap(jet);
hp4 = get(subplot(2,2,4),'Position');
if max(max(max(abs(funcM))))>0.1
    colorbar('Ticks',round(min(min(min(funcM))):(max(max(max(funcM)))-min(min(min(funcM))))/5:max(max(max(funcM))),2),...
             'Position', [hp4(1)+hp4(3)+0.01  hp4(2)  0.03  hp4(2)+hp4(3)*2.1])
else
    colorbar('Position', [hp4(1)+hp4(3)+0.01  hp4(2)  0.03  hp4(2)+hp4(3)*2.1])
end

ax=axes('Parent',f,'visible','off'); 
xlim(ax,[0,1]);ylim(ax,[0,1]);
xticks(ax,0:0.2:1);yticks(ax,0:0.2:1);
ax.Title.Visible='on';
ax.XLabel.Visible='on';
ax.YLabel.Visible='on';
ylabel(ax,'Position X/L');
xlabel(ax,'Position Y/L');
title(ax,['OD = ',num2str(obj.physPara.OD),' ,P_{in}/P_* = ', num2str(obj.physPara.power/obj.physPara.P0)])

end

